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Part I. I i ine E’i€..ien^ 

A.' Time Transformations 

Time transformations of the Sundman type (ref. 1) have the form 

dt = cr a ds (1 ) 

where t Is the time, c and a are positive constants, r the magnitude 
of the radius-vector, and s the new independent variable. 

It has been known for many years that use of the time transformations 
of equation (1) improves rates of convergence of analytical series solutions 
of gravitational systems (ref. 2 and 3). Also, when equation (1) is used 
in conjunction with a coordinate transformation, singularities due to col- 
lisions may be eliminated from the equations of motion (see, for example, ref 4) 

Several recent published studies (ref. 5, 6, 7) and unpublished studies 
(ref. 8, 9) show that the use of the time transformations given by equation (1) 
with 1 < a < 2, substantially improves numerical integration accuracy and 
efficiency for satellite solutions. Moreover, it has been shown by many of 
these same studies that, for a large class of satellite initial conditions, 
a near or equal to 1.5 provides the best accuracy and efficiency (ref. 7,8,9). 

B. Time Elements 

It appears that time transformations lessen or weaken to some extent 
the dynamical in-track (Liapunov) instability associated with satellite motion. 
The reduction of instability in the coordinates is rendered at the expense of 
the introduction of a differential equation in order to obtain the time 
(equation (1)). The differential equation for the time exhibits some of the 
dynamical instability that was removed from the differential equations for the 
coordinates. If the perturbations are small relative to the two-body forces. 
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then souse the instability of the time equations may be removed by use of 
time elements. 

With zero perturbations, equation (1) may be integrated in closed form 
for o = 1, 1.5, 2. These integrals define "time elements" (ref. 11). For 
perturbed motion, the integrals of equation (1) are then differentiated using 
methods of variation of parameters. 

The resulting differential equations for the time elements are then 
integrated by means of numerical integration. A significant portion of the 
instability is removed from the equation that defines the time transformation 
(equation (1)) by taking account of the two-body solution. 

It has been shown by Stiefel and Scheifele (ref. 11) that a time element 
may be introduced for the KS equations of motion in the four dimensional u-space 
The related time transformation has a - 1 in equation (1). 

Often in practice, one does not wish to use the complex KS coordinate 
transformations to four dimensional space, but rather remain the three 
dimensional physical space and use only the time transformation of equation (1). 
Also, the time transformation for optimum accuracy has a * 1,5 in equation (1) 
For these reasons we have derived a time element for a = 1.5, as well as for 
a = 1 and 2, in the three dimensional physical space. These time elements and 
their defining differential equations are presented here in this report. 


C. The New Time Elements for a « 1, 2 
For two-body motion and for a = 1 and 2, the integral of equation (1) 
is Kepler's equation for certain choices of the constant c. 

For a = 1 and c = /~~ , the integral is 


t = 


T 



e sin E, 


( 2 ) 


where E is the eccentric anomaly, and where a and e are the semi major 



axis and eccentricity, respectively, n.nd t is the time element. A 
differential equation for the time element is 
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• d_r 
ds 



, 1 . - de_ 3 . .i„ P A 1 da 

r /y sin E ds 2 sin E /y ds 

+ a 3 e cos E { (1 -e cos E)vE - sin E ve} • R 


(3) 


where vE, and ve denote the gradients of E and e with respect to the 
velocity vector, and R is the perturbing force. 

The first term on the right hand side of equation (3) is not a constant 
for perturbed motion and is of the - order of magnitude of the unperturbed 
quantities. The semi -major axis in the term must be determined from the state 
vector which has errors due to the numerical integration of the equations of 
motion. The subsequent integration of these errors in equation (3) to deter- 
mine the time element compounds the error substantially. To reduce this large 
source of errors in the time element the following relation may be used to 
replace the first term on the right in equation 3: 



_ _P - 

/ -8(h-v)3 ' 


(4) 


where h is the total energy of the perturbed system aid V is the perturb- 
ing potential defined by 


h - h k + V 

where is the total unperturbed or Keplerian energy. 
For a - 2 in Equation (1), and 

c = — L_ 

/ya(l-e z / 


the integral is the same as equation (2) but here s is equal to the true 
anomaly, f. The differential equation for the time element is 
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dx _ ar 
ds v£p 



a(l-e 2 ) 




sin E ds 


rr** 
/a 5 e 

/ \i r 


sin E cos E g|- + a 2 re cosE(l-e 2 )~^ 


*{(l-e cos E) vE - sin Eve} • R (5) 

A substitution similar to that of equation (4) may be used to replace part of 
the first term on the right in equation (5). 


D. The Intermediate Anomaly 

For a = 1.5 and c = 1, the integral of equation (1) for two-body 
motion is 


s 


2 

Trmry 


F{f, k) 


( 6 ) 


where s is the new independent variable, referred to here as the "intermediate 
anomaly," e is the eccentricity, F the incomplete elliptic integral of the 
first kind, k the modulus defined in terms of the eccentricity as 



and f the true anomaly. Equation (6) may also be written as 

sin -^f = sn[|- /1+e s] (7) 

where sn is the sine Jacobian elliptic function. 

Equations (6) and (7) give the true anomaly as a function of s. The 
time is given by equation (2) where the eccentric anomaly is obtained from the 
true anomaly. The time element is obtained from the following differential 
equation. 
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dx 
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ar 1 /^ + ~ / — e sin u ~r + / — 
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^ de_ 
ds 


rsm u - e cos u 


tan — f (1 + cos u) 


( l + e)/!^ 2 


-i-/£ 

2e a 
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2E 


(1+e) 1 / 2 (l-e) 


Sn(^s/Tl+eTp") cn (~ $/(ire)p 
dnC—S/l+eJy ^ 


( 8 ) 


The quantity E Is the incomplete elliptic integral ef the second kind; and 
sn, cn, and dn are the sine, cosine, and amplitude Jacobian elliptic 
functions, respectively. 

The differential equations for any of the time elements that we have 
derived have no pure or mixed secular terms for unperturbed motion, therby 
further reducing the in- track (Liapunov) instability associated with 
equation (1) 


E. Elliptic Function and Integral Algorithms 
We have developed new computational algorithms tailored to our purpose 
in order to compute the necessary complete and incomplete elliptic integrals 
and the Jacobian elliptic functions. The algorithms are based on and 
modifications of those by Hofsommer and van de Riet (ref. 12) and are about 
five times faster than the algorithms of DiDonati and Hershey (ref 13). 


F. Conclusion 

We have compared the accuracy of integration of equation (1) directly with 
that of integrating equation (3) using equation (2). The result is that use 
of equations (2) and (3) provide about two to three times the accuracy of equa- 
tion (1) after about three satellite revolutions. Further studies are presently 
being performed for longer integration times and for equations (5) and (8). 

The full paper with a complete set of numerical comparisons will be pre- 
sented at the COSPAR Symposium on satellite dynamics, June 19-21, 1974, 



Sao Paulo, Brazil, and at the AIAA/AAS Astrodynamics Conference, August 5-9, 
1974, in Anaheim, California. The paper will be submitted for publication to 
the Journal of Celestial Mechanics or be published in the Proceedings of 
the COSPAR meeting. 

Part II. Stabilization by External Energy Corrections 

A. INTRODUCTION 

It is shown in reference (14) that constraining solutions to satisfy 
exactly the integrals of motion tends to Liapunov stabilize the solutions of 
unstable dynamical systems. It is shown in reference(15) , pp. 73-77, 
that complete Liapunov stabilization of Keplerian motion may be realized by 
regularization. Reference (15) indicates that the regularized equations of 
motion are stable due to the fact that the Keplerian energy is constant and 
that the constant appears explicitly in the differential equations. This 
is complete stability in the sense of Liapunov since the time equation is also 
stabilized along with the dependent variable equations. This point confirms 
the conclusions of reference ,(14) since reference (14) shows that use of 
integrals in stable dynamical systems (harmonic oscillator) does not improve 
accuracy. Whereas use of integrals for unstable (in the sense of Liapunov) 
dynamical systems (for example, Keplerian motion) improves the accuracy of 
the solution by several orders of magnitude. In addition, references (14) 
and (16) indicate that neither time nor coordinate transformations are nec- 
essary for the stabilization, only the use of the integrals of motion. 
Reference (15), p. 76, reference (16), reference (17), and reference (18) 
indicate that for unperturbed and perturbed Keplerian motion, only the energy 
integral is needed for stabilization. This result is in agreement with the 
numerical results of reference (14), which show that for an n-body problem 
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dynamical system, by satisfying only the energy integral exactly, a solution 
is produced that is several orders of magnitude more accurate than by not 
satisfying the integral and very nearly as accurate as a solution satisfying 
all ten integrals of motion. Hence, it appears that there is a fundamental 
relation between isoenergetic solutions, regularization, and dynamical 

V 

(l^apunov) stabilization. 

In reference (14), corrections are applied to the components of the 
state vector so that the corrected state satisfies identically the integrals 
of motion. The corrections are chosen so that the sum of the squares of the 
corrections is a minimum. 

This report applies the concept of stabilization using integrals to the 
solution of the equations of motion of artifical satellites in an attempt to 
reduce the propagation of local truncation errors and improve the efficiency of 
a numerical integration solution process. Even though no constant energy 
integral exists for the full satellite equations of motion due to the presence 
of the tesseral harmonics, drag, and third-body perturbations; nevertheless, 
a slowly varying energy "integral" is present in a coordinate system rotating 
with the rotating Earth. The integral in the rotating coordinate system is 
analogous to the Jacobian integral in the restricted problem of three bodies. 

The method presented in reference (14) is extended here to the satellite 
solution in a rotating coordinate system. An application of the method to 
satellite solution is presented showing accuracy improvements of at least 
two orders of magnitude with negligible additional computation time. In 
addition, a modification of the method presented in reference (14) is presented 
to allow the use of slowly varying integrals. 

B. The Equations of Correction 

The equations of correction using integrals derived in reference (14), 
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will now be adapted to a dynamical system with three degrees of freedom 
using only one integral - the energy intergal. 

For a dynamical system with three degrees of freedom, let x = [x-j ,x 2 ,x 3 »x^ ,Xg,Xg] 
be the state vector in the phase space, where and x^ are the 

coordinates and x^, x^ and Xg are the corresponding velocity components. 

Let 

E{x) = 0, (1) 

be an integral of the system. Equation (1) defines a hypersurface of five 
dimensions imbedded in the phase space of six dimensions. 

During a process of numerical integration of the system, a computed 
solution is obtained at time t: 

n ~ n(t) *Tig >Tigl $ 

where n-| » and n 3 are the computed position components and and 

Tig are the computed velocity components. Due to errors in the computational 
procedure, the integral may not be satisfied exactly but 

E{n) = e » (2) 

where e is the small quantity. The solution has left the integral surface 
defined by Equation (1) and is on the surface defined by Equation (2). It 
is desired to make corrections An = [at^ »An 2 ,An 3 ,An 4 ,An 55 Ang] to the 
computed vector n to obtain the vector 

x = n + An , 

such that 

E(x) = 0 . (3) 

The square of the magnitude of the correction vector An may be written as 
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f(fin) = I (fin,) 2 . (4) 

1=1 

The corrections are uniquely chosen so that the function f of Equation (4) 
is minimized, subject to the constraint of Equation (3). The solution may be 
obtained by solving the following equation 


An i - = 0 , i = 1,2,3, 4,5,6 (5) 

along with Equation(3) for the seven unknowns A, and An^- » i * 1,2,3, 4, 5, and 
6. The quantity A is the Lagrange multiplier. Unless the integral given 
by Equation {3} is a simple function of the variables (for instance linear), 
the solution of the system may be complex (or. perhaps not obtainable) as is 
the case when the integral is the integral of energy of a gravitational 
system. The solution may be simplified by an expansion of the integral in 
powers of the corrections. The expansion becomes 

E(x) = E(n) + I — in. + ... . (6) 

1=1 an 1 

Since the errors of the computation and hence the necessary corrections, An^ > 
are generally small, second and higher-order terms may be neglected. 


Solving Equations (5) and (6) for the correction An., with E(x) = 0 
and E(n) = e» yields ;* . 




f T3E 


2 


t — 1,2, 3, 4, 5, 6. 


(7) 


The correction vector An is added to the computed state vector n to 
obtain a new state vector x which satisfies the integral E(x) = 0, with 
an error of order j An | . Geometrically, minimizing Equation (4) subject 
to Equation (6) causes the vector An of Equation (7) to be normal to a 
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five-dimensional plane which is approximately tangent to the surface 
E = 0 at the point x. The equation of the plane is given by Equation (6), 
neglecting the second and higher-order terms. Equation (7) is generalized 
in reference (1) to a system having n degrees of freedom and p integrals 
of motion, where n and p are any positive integers. Also, a more detailed 
derivation of Equation (7) is given in reference (14). 

The energy integral of a dynamical system (per unit mass) may be written 
in terms of a potential energy U and the state vector x as 

E(x) - \ (xj + + x*) + U - C = 0 

where C is the value of the energy for a set of initial conditions. 

A numerical integration of the system yields the solution vector n 
at time t. The errors in the computation may cause E to be nonzero: 

E(n) - £> (8) 

where e is the error of the integral. The correction vector An may be 

calculated by using Equation (7), so that E(n + An) = 0. For the calculation, 

3E 

the quantities e and - — are needed. 

3n_- 

3E 1 

The quantities - — are the partial derivatives of E with respect to 
n(or x) and are easily computed. The partial derivatives of E-j with 
respect to x^, Xg, and Xg are equal to x^, Xg, or Xg, respectively. 

The partial derivatives of E-j with respect to x-j , X 2 » and are 
simply minus the respective components of the force F, already computed 
during the integration step. 

C. Previous Applications and Evaluations of the Method 
As presented in reference (14), the method was applied to the numerical 
integration of several dynamical systems to determine its paractical value. 

The systems considered were the harmonic oscillator, the gravitational system 
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of two-bodies, and the gravitational system of 25-bodies. 

The application of the correction method to the harmonic oscillator in 
a phase space of two dimensions showed no noticeable differences in accuracy 
between the corrected and uncorrected solutions (reference 1). 

However, the application of the method to the system of two-bodies 
in a phase space of four dimensions over a range of initial conditions showed 
a large difference in accuracy between corrected and uncorrected solutions. 

The corrected solutions were about three orders of magnitude more accurate than 
the uncorrected solutions (reference (14)). 

The different results obtained for the harmonic oscillator and the system 
of two-bodies offers an explanation of when and why the method appears to be 
of value, as presented in reference (14). The errors in the integral of the 
harmonic oscillator are small compared to the error in the state variables 
of the solution. Since the harmonic oscillator is a stable system, a solution 
with a small error will not diverge from a system without the error. This 
would explain the result indicating no difference between the corrected and 
the uncorrected solutions for the harmonic oscillator. The errors in the 
energy integral of the system of two-bodies are also small relative to the 
total errors in the state variables of the solution. But the system of 
two-bodies is unstable in the Lyapunov sense and hence the system with the 
energy errors will diverge from the system without the energy errors. 

The method was applied to a gravitational system of 25-bodies. The 
results show that the method yields a more efficient numerical integration 
process of the n-body problem. A greater accuracy of about two orders of 
magnitude is obtained with the method while using the same time of calculation 
then without using the method. Or the same accuracy is obtained with the 
method while using about 25% less time of calculation (reference 14). 



12 


0. App lication of the Method to Satellite Motion 
Recently, the method has been applied to the solution of an artificial 
satellite of the Earth. The equations of motion of the satellite included all 
zonal, sectoral, and tesseral spherical harmonics. The equations of motion 
did not include luni-solar nor drag effects. With all of the spherical 
harmonics included, no energy integral exists. However, an integral exists 
in a coordinate system rotating with the rotation of the Earth. The integral 
is analogous to the Jacobian integral in the rotating coordinate system of the 
restricted problem of three bodies. If the integral is formulated in the 
rotating coordinate system using position and velocity relative to the rotat- 
ing coordinate system and then transformed into the fixed system, in terms of 
position and a velocity relative to the fixed system, the integral has the form 

E(y y) = ^ y*Y +U-C+y®Y* W - 0 (9) 

where U is the potential, y> is the position vector, and y is the velocity 
vector relative to the fixed frame. The quantity C is the value of the 
integral for a certain set of initial conditions and W is the angular velo- 
city vector of the rotation of the Earth. The quantity y & y • W is the 
component of the angular momentum in the direction of the rotational axis of 
the Earth times the rotation rate of the Earth. That equation (9) is a 
constant of the motion may be proved in several ways. One way is to formulate 
the Hamiltonian of the satellite system in extended phase space. Equation (9) 
is that Hamiltonian, neglecting the constant C. 

In addition, a modification of the method presented in reference (1) may be 
be incorporated to allow the use of slowly varying integrals. A differential 
equation may be formulated giving the rate of change of an integral. This 
equation may be integrated numerically along with the state vector. Then, 
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at any time, the value of the integral is known and the correction procedure 
of reference (1) is applied. It appears that such a modification induces 
stabilization and does so only for slowly varying integrals. 

E. Numerical Results 

f 

The method was applied to the satellite solution by the following procedure. 
Two solutions were obtained by numerical integration. One solution did not 
utilize the integral while the other introduced corrections determined by 
equations jj ) , (8), and (9). The corrections were applied to the state vector 
after each integration step, as in~the discussion following Equation (1) 
above. Both solutions were compared with- a more accurate solution determined 
by smaller integration step-sizes. Both solutions were obtained using a high- 
order predictor-corrector integration algorithm of the STDS system of the NASA 
Goddard Space Flight Center (reference 19). Both uncorrected and corrected 
solutions were integrated with the same step-size and with the same number of 
integration steps. 

The satellite had initial orbital parameters as follows: A semimajor axis 
of 136000 km; an orbital period of 120 hours; an eccentricity of .8; and an 
inclination of 40°. The solutions were found for about 5 revolutions or over 
a time span of 3 weeks. The corrected solutions were found to be more accurate 
than the uncorrected by at least two orders of magnitude. That is, the position 
and velocity components of the corrected solution had at least two more digits of 
agreement with the true solution than the uncorrected solution. 

In addition, solutions obtained by a 4th-order Runge-Kutta were performed 
for a satellite with an eccentricity equal to 0.6 and a period equal to 3 hours. 
Two-body forces and forces due to the second zonal harmonic were included. 
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Various step sizes were used: from 500 steps per revolution to 4000 

steps per revolution. After 10 revolutions, the corrected solution was at least 
one, and often two, orders of magnitude more accurate than the uncorrected 
solution,- for all step sixes used. 

Also, the following study was performed. Corrections were made according 
to the following two techniques: (1) corrections were made to both the 

position and velocity components with various step sizes; and (2) corrections 
were made to only the position components. This study was initiated so as: 

0) to study whether there is a difference if corrections are used with a 
Class I integrator or with a Class II integrator; and (2) to determine whether 
reapportioning the corrections from the velocity vector to the position vector 
could cause differences. 

Result: No significant difference was found between the corrected final 

state for technique (1) and for technique (2). 

Part III. Long-term Global Solutions for the Synchronous Satellite. 

A preliminary report. 

A. Introduction 

A new method has been developed by one of the principal investigators 
(Nacozy) to yield semi -analytical solutions of the long-term solution of reso- 
nant satellites, in deep resonance with the tesseral harmonics. 

The method is able to yield global solutions for all eccentricities, 
inclinations, and commensurability ratios for resonant satellites, including 
the synchronous satellite. The method is not restricted to moderate or small 
eccentricities and inclinations nor to just synchronous satellites as is the 
analytical solution of Musen and Baile (reference (21)). In addition, the 
method does not yield particular solutions as are given by Gedeon (reference 
(22)) and others. 
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Our method is based on a semi -analytical approach to resonant pertur- 
bations introduced by Hori (reference (23)), developed by Giacaglia (reference 
(24)), and extended and applied by Giacaglia and Nacozy (reference (25)), 
and Nacozy and Diehl (reference (26) (and partially used by Musen and Bailie 
reference (21)). The method numerically averages the Hamiltonian and uses 
the fact that the averaged Hamiltonian has a minimum with respect to the 
resonant (critical) argument at the stable, stationary value of the resonant 
argument. 

Our semi-analytical method also appears to have the advantage that luni- 
solar and radiation pressure effects may be easily added numerically without the 
need for analytical developments. 

We are presently applying the method to synchronous satellites, to determine 
the practicality and potential of the method. For the ^ 2~^22 we have 

now successfully obtained long-term solutions for e, i, u, and n, 
for a]j_ e and i, for a synchronous satellite in the vicinity of the stable 
stationary solution. We have discovered a family of periodic solutions for 
the ^ 2~^22 P ro ^ em * These orbits are inclined and eccentric orbits and 
have a repeating ground track. 

We next plan to include more zonals and Tesserals and luni-solar perturbations 
and radiation pressure. We also have the capability to consider the 4, 6, 12, 
and IS hour resonant satellites. 

Part IV. The D-S Equations of Motion 

Abstract 

A new set of canonical two-body elements referred to as Delaunay-similar 
(D-S) elements is presented in references (27) and (28). In contrast to the 
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classical Delaunay theory which has time as the independent variable, the D-S 
theory uses an independent variable which is a generalized true anomaly. This 
study^is concerned with numerical integration of the canonical perturbation 
equations of these elements. A description of the derivation of the D-S 
elements is given. Two modifications are introduced which increase the 
numerical stability of the system. The differential equations are established 
in Gaussian form fit for numerical integration. All associated transformation 
formulae and partial derivatives are described in detail. 

(This is an abstract of a study performed by G. Scheifele and R. Samway. The 
complete study is given in the Masters Thesis of Robert C. Samway, The 
University of Texas at Austin, August, 1973. The Thesis may be obtained 
through the Department of Aerospace Engineering and Engineering Mechanics, 

The University of Texas at Austin.) 
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